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Abstract. We present calculations of the spectral and temporal radiative signatures expected from ultrarelativistic protons in 
compact sources. The coupling between the protons and the leptonic component is assumed to occur via Bethe-Heitler pair pro- 
duction. This process is treated by modeling the results of Monte-Carlo simulations and incorporating them in a time-dependent 
kinetic equation, that we subsequently solve numerically. Thus, the present work is, in many respects, an extension of the lep- 
tonic 'one-zone' models to include hadrons. Several examples of astrophysical importance are presented, such as the signature 
resulting from the cooling of relativistic protons on an external black-body field and that of their cooling in the presence of 
radiation from injected electrons. We also investigate and refine the threshold conditions for the 'Pair Production/Synchrotron' 
feedback loop which operates when relativistic protons cool efficiently on the synchrotron radiation of the internally produced 
Bethe-Heitler pairs. We demonstrate that an additional component of injected electrons lowers the threshold for this instability. 
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1. Introduction 

The spectral energy distribution (SED) of powerful AGN such 
as flat-spectrum radio quasars and blazars has a double humped 
appearance with the low energy part extending from the radio 
to UV (or in extreme cases to X-rays), and a high energy part 
extending from X-rays to y-rays. In AGN with relativistic jets 
closely aligned to the line of sight the emission is dominated by 
non-thermal radiation, with the low energy hump being mainly 
synchrotron radiation. If the alignment is not so close, a thermal 
component of UV radiation from an accretion disk may dom- 
inate. The non-thermal components can be strongly variable, 
probably originating in the jet. 

These observations indicate that the jets of blazars act as 
eflicient particle accelerators. Furthermore the gamma-ray ob- 
servations in the GeV (Hartman et al. 1999) and TeV regime 
(Horan & Weekes 2004) imply that the accelerated particles 
can reach very high energies. Models involving electron radi- 
ation can adequately explain both this high energy emission 
and the coordinated multiwavelength campaigns (Mastichiadis 
& Kirk 1997; Tavecchio et al. 2001; Krawczynski et al. 2002). 
The usual assumption is then that the high energy part of the 
SED is due to inverse Compton scattering of the low energy 
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part of the SED (synchrotron self-Compton model) possibly 
supplemented by inverse Compton scattering of external pho- 
tons (external Compton models) for example from the disk, ei- 
ther directly or scattered by clouds. Despite these successes, 
the question of the role of a possible relativistic hadronic com- 
ponent remains an open one. 

In principle, sites of electron acceleration may accelerate 
protons as well. Consequently, models in which the high energy 
part, and some fraction of the low energy part, of the SED is due 
to acceleration and interaction of protons in the jet have also 
been proposed. Some of these models invoke interactions with 
ambient matter (Beall & Bednarek 1999; Pohl & Schlickeiser 
2000; Schuster et al. 2002) but they require high mass densities 
in the jet to be viable. Here we concentrate on hadronic mod- 
els in which the protons interact with low energy photons via 
Bethe-Heitler pair production. 

As with leptonic models, the target photons may be pro- 
duced inside the emission region in the jet or may originate 
from outside of the jet, e.g., from an accretion disk (Protheroe 
1997; Bednarek & Protheroe 1999; Atoyan & Dermer 2001; 
Neronov & Semikoz 2002). For internally produced target 
photons, synchrotron emission by a co-accelerated popula- 
tion of electrons is assumed (Mannheim & Biermann 1992; 
Mannheim 1993, 1995). The high energy hump of the SED 
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then results from electromagnetic cascading of gamma-rays 
from tt" decay and electrons from Bethe-Heitler pair produc- 
tion and ;r'^ ^ ju^ ^ e"^ decay in the radiation and mag- 
netic field of the blob. Neutrinos and cosmic rays would also 
be emitted as a result of the neutrinos from the n- and jj.^ 
decays, and neutrons produced in pj — > nn"^ interactions if 
the threshold for pion photoproduction is exceeded (Eichler 
& Wiita 1978; Sikora et al. 1987, 1989; Kirk & Mastichiadis 
1989; Begeknan et al. 1990; Giovanoni & Kazanas 1990; 
Protheroe & Szabo 1992; Szabo & Protheroe 1994; Waxman 
& Bahcall 1999; Mannheim et al. 2001; Atoyan & Dermer 
2003; Protheroe 2004). For protons to be accelerated to en- 
ergies sufficient to exceed the Bethe-Heitler and photo-pion- 
production thresholds, relatively high magnetic fields are re- 
quired. In proton synchrotron blazar models the magnetic field 
is sufficiently high such that the the high energy part of the 
SED has a major contribution also due to synchrotron radiation 
by protons (Mticke & Protheroe 2000, 2001; Aharonian 2000; 
Reimer et al. 2004). All of the work described above assumes 
the emission has reached a steady state, and that the target pho- 
ton fields are steady. In reahty, the strong variabihty displayed 
by these sources mandates a time-dependent calculation, that, 
ideally, should be done self-consistently, with internally pro- 
duced radiation fields contributing alongside external ones to 
the target radiation field. 

Time-dependent codes that solve the kinetic equations de- 
scribing electrons and photons and their interactions have been 
developed and successfully appUed to AGN (Mastichiadis & 
Kirk 1997; Krawczynski et al. 2002). However, codes of this 
type that also account for hadronic interactions have been ne- 
glected. One reason for this is that whereas the modeling of 
leptonic processes is relatively straightforward (e.g., Lightman 
& Zdziarski 1987; Coppi & Blandford 1990), photo-hadronic 
and hadron-hadron interactions are much more complex. To 
date, all attempts have used approximations of uncertain accu- 
racy (e.g.. Stern & Svensson 1991) but the use of Monte-Carlo 
event generators which model in detail electromagnetic (Szabo 
& Protheroe 1994; Protheroe & Johnson 1996) and hadronic 
interactions (Miicke et al. 2000) opens up the possibility of ex- 
tracting accurate descriptions of the fundamental interactions 
suitable for incorporation into a kinetic code. 

Motivated by these developments we investigate the con- 
sequences of the presence of relativistic hadrons in compact 
sources by incorporating new results from Monte-Carlo simu- 
lations into a time-dependent code which follows the evolution 
of relativistic hadrons, electrons and photons by solving the 
appropriate kinetic equations. In the present paper we inves- 
tigate as a first step, the case in which the only channel of cou- 
pUng between hadrons and leptons is the Bethe-Heitler pair- 
production process, leaving the investigation of photo-meson 
production for a future paper. Although this is not a complete 
description of hadronic models it nevertheless enables one to 
draw useful and interesting new results. 

The present paper is structured as follows: In Section 2 we 
present the numerical code that solves simultaneously in a self- 
consistent manner the coupled, time-dependent kinetic equa- 
tions for each species, i.e. protons, electrons and photons. In 
Section 3 we present the Monte-Carlo results for the Bethe- 



Heitler process and show how these can be incorporated in 
the kinetic equations. In Section 4 we present some results for 
the case in which relativistic protons interact with an external 
black-body radiation field. In Section 5, we present a numer- 
ical analysis of the 'Pair-Production/Synchrotron' instability. 
The case of simultaneous injection of relativistic protons and 
electrons is examined in Section 6 and the main conclusions 
are summarized in Section 7. 

2. The kinetic equations for electrons, protons and 
photons 

The kinetic equations describing a homogeneous source region 
containing protons, electrons and photons were formulated and 
solved numerically by Mastichiadis & Kirk (1995, henceforth 
MK95). We follow the same method, using an improved de- 
scription of the microscopic processes. The equations to be 
solved can be written in the generic form 

-^+Li + Qi = (1) 
ot 

where the index i can be any one of the subscripts 'p', 'e' or 
'y' referring to protons, electrons or photons respectively. The 
operators L,- denote losses and escape from the system while Qi 
denote injection and source terms. These are defined below. 

The unknown functions «, are the dififerential number den- 
sities of the three species, normalised as follows: 

Protons: 

n*p(yp,t)dyp = (rTRnp(Ep,t)dEp withyp = (2) 

Electrons: 

Et 

nl(ye,t)dye - <TTRne(Et,t)dEe with ye = r (3) 

Photons: 

n*,{x,t)dx = (T'YRny{ey,t)dey withx = ^ (4) 

and the time t has been normalised in all equations to the light- 
crossing time of the source ta = R/c. 

The physical processes to be included in the kinetic equa- 
tions are: 

1 . Proton-photon (Bethe-Heitler) pair production which acts 
as a loss term for the protons (L^^) and an injection term 
for the electrons (Q^^) 

2. Synchrotron radiation which acts as an energy loss term for 
electrons (Lf^) and as a source term for photons (2y^) 

3. Synchrotron self absorption which acts as an absorption 
term for photons (L^'^) 

4. Inverse Compton scattering (in both the Thomson and 
Klein-Nishina regimes) which acts as an energy loss term 
for electrons iV^^) and as a source term for high energy 
photons and a loss term for low energy photons, both ef- 
fects included in Q^^^ 

5. Photon-photon pair production which acts as an injection 
term for electrons (Q^) and as an absorption term for pho- 
tons (L^) 

6. Electron-positron annihilation which acts as a sink term for 
electrons (Lf^) and as a source term for photons (Q^) 
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7. Compton scattering of radiation on the cool pairs, which 
impede the free escape of photons from the system. This 
effect is treated approximately by multiplying the photon 
escape term by the factor (1 +Hil -x)tt/3)"\ where tj the 
Thomson optical depth, while H(l-x) is the step-function. 
(Lightman & Zdziarski 1987) 

With the inclusion of the above terms the kinetic equations 
become, for each species (from now on we refer only to nor- 
malised quantities and, for convenience, drop the asterisks) 

- Protons 

t/t *p,esc 

- Electrons 

^ + Lr + 1}:' +Lr+^ = qt + + er (6) 

OT 'e,esc 

- Photons 
fin fi 

any ny ^ . 

dt 1 -I- Tt/3 ^ y ^^y ^^y W 

We note the following regarding the above equations 

1. When the various terms above are written explicitly, equa- 
tions (5), (6) and (7) form a non-linear system of coupled 
integro-differential equations. 

2. The various rates conserve the energy exchange between 
the species - for example, the amount of energy lost per 
second by electrons at each instant due to synchrotron radi- 
ation is equal to the power radiated in synchrotron photons, 
(for details see MK95) 

3. In the absence of the Bethe-Heitler pair-production term 
2™, equations (6) and (7) decouple from the protons 
(Eq. 5) and the system becomes identical to the 'one-zone' 
time-dependent leptonic models 

4. Protons are injected via the term with a prescribed dis- 
tribution in energy. Thus, in contrast to MK95, we do not 
investigate the effects of particle acceleration. 

5. Electrons may also be injected externally through the pre- 
scribed term ^e"'- However the two other injection terms in 
the electron equation {(^^ and (^J) are determined self- 
consistently from the proton and photon distributions. 

6. Both protons and electrons (we make no distinction be- 
tween electrons and positrons in the present treatment) es- 
cape from the source region on the timescales ?p,esc and ?e,esc 
respectively (given in units of t„) 

We can also define various compactnesses related to the 
most important of the above quantities. So we define the photon 
compactness as (see MK95) 

1 r ny(x, t) 
(y = - \ dxx (8) 

while the scaled to electron rest-mass compactness of exter- 
nally injected protons is 

= ^ r dy{y-\)Q^{y). (9) 
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In analogous fashion we can also define the compactness of the 
externally injected electrons 

C - \ ^ djij-DQTiy). (10) 

Finally in order to calculate the Thomson optical depth of the 
cool pairs we use 

X1.26 
dyrieiy), (11) 

Closing this section we note that the treatment of syn- 
chrotron and inverse Compton scattering has been improved 
over that described by MK95 in that the full emissivities 
are incorporated, rather than delta-function approximations. 
However, the main improvement is in the treatment of the 
Bethe-Heitler pair-production process, using Monte-Carlo sim- 
ulations, as described in the following section. 

3. Bethe-Heitler pair production 

3.1. Monte-Carlo simulations 

For an isotropic target comprising monoenergetic photons of 
energy e = xrricC^ the effective cross-section for interaction of 
a proton of energy E = y^nipC^ is given by 

(o-BH(rp' •^)> = 2 J -/3p COS e)cr(s)d COS (12) 

where is the angle between the proton and photon directions, 
^min(7'p, x) is the minimum value of this angle consistent with 
the threshold, 

s = m^c'* + 2e£'( 1 - ySp cos 0) (13) 

is the centre of momentom (CM) frame energy squared, ySpC 
is the proton's velocity, and ctbh the total cross-section for 
which we use the Racah formula as parameterized by Maximon 
(1968) (see Formula 3D-0000 in Motz et al. 1969). Changing 
variables, one obtains the angle-averaged cross-section 

(o-BH(yp,x)) = I (Tis)is -mlc'^yds, (14) 

where 

= (mpc^ + InieC^f « 0.882 GeV^, (15) 

and 

*max(rp> X) = mlc'^ + 2ypOTpC'XmeC^(l +/?p) (16) 

corresponding to a head-on collision. The angle-averaged 
cross-section is plotted as a function of the product of photon 
energy and proton energy in Fig. 1. For x <c 1, the threshold 
condition impUes 7p » 1 and j8p » 1, so that to a first approx- 
imation, the cross-section is a function of xjp, rather than of x 
and 7p separately. 

Examination of the integrand in Equation 14 shows that the 
square of the total CM frame energy is distributed as 

p{s) oc cr(i)(i - m.pc'^), (17) 
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Fig. 1. The angle-averaged 
section (crBnCyp,^)) (solid 



Bethe-Heitler pair-production cross 
curve), the angle-averaged Bethe- 



Heitler pair-production inelasticity ^bhCTp,-*) (dashed curve), and 
^bh(Tp, xXa-Bniyp, x)) (dotted curve) are plotted as a function of y^x 
for X <s 1. 



in the range .Vmin < v < Vmax- The Monte Carlo rejection 
technique is used to sample s, and Equation 13 is used to 
find 6. We then Lorentz transform the interacting particles to 
the proton rest frame and sample the positron's energy from 
the single-differential cross-section, do-/dE+, for which we 
use the Bethe-Heitler formula for an unscreened point nucleus 
(Formula 3D-1000 in Motz et al. 1969). Finally, the positron's 
direction is sampled from the double-differential cross-section, 
d(r/dE+dCl+ for which we use the Sauter-Gluckstem-HuU for- 
mula for an unscreened point nucleus (Formula 3D-2000 in 
Motz et al. 1969), and its laboratory frame energy is obtained 
by a Lorentz transformation. For a range of proton energies, 
the simulation is repeated a large number of times to build up 
distributions in energy of positrons produced in BH pair pro- 
duction. 

The distribution in energy y^m^c^ of electrons (of either 
charge). 



/(re;rp,x)-^, 



(18) 



is taken to be twice that for positrons, as discussed in Protheroe 
& Johnson (1996), and is plotted in Fig. 2 for x = 10"^ and 
three yp values. In a similar calculation using black-body target 
photons (Protheroe & Johnson 1996), the mean inelasticity was 
found to be in excellent agreement with those calculated ana- 
lytically (Blumenthal 1970; Chodorowski et al. 1992; Rachen 
& Biermann 1993). The mean inelasticity, given by 



^(7p,x)^— I dye—f(ye;yp,x) 
mp J yp 



(19) 



is plotted as a function of the product of photon energy and 
proton energy in Fig. 1 (dashed line) together with its product 
(dotted curve) with the cross-section, to show how the proton 
energy loss rate depends on energy. 



Fig. 2. The distribution of electron energies is shown for x= 10"* and 
yp = 10*-^ (dashed histogram), 10^ (soUd histogram) and 10' (dotted 
histogram). 



3.2. Numerical Modeling 

3.2.1. Electron/positron Pair Production 

The pair-production spectra were calculated using the results 
of the Monte Carlo code described in Section 2.1. Protons of a 
specific energy were allowed to interact with isotropic monoen- 
ergetic target photons and the energies of the products were 
tabulated. The photon target energies used were xq = 10"*, 
10"^, 10"^ and 1 (in units of electron rest mass). Proton ener- 
gies ranged from y^ = x^^ (so the threshold requirement could 
be met) up to yp = 10"^ Xq' in logarithmic steps of 0.1. 
The pair-creation rate is then given 

2e"(re,0 = J dy^HpiypJ) J dxnyix,t)x 

f(yp,yc,x){(rBHiyp,x)) (20) 

where /(7p,7e. ^) is the distribution found from the Monte- 
Carlo modelling (see eq. 18), normalised such that 



dyj{yp,ye,x) = 2. 



(21) 



As simple parameterisations of these curves do not give ac- 
ceptable fits, we tabulated the spectra and used interpolations 
to derive the produced pair-injection spectrum from a particular 
proton-photon collision. 

We note that MK95 used the approximation fiyp,ye,x) - 
2<5(7e - Tp) corresponding to f = /We/mp. For a proton distribu- 
tion np(yp) which decreases as yp increases, this overestimates 
the production rate of electrons of a given energy. 

3.2.2. Proton losses 

Since the proton loses a small amount of energy (typically 

given by Ayp ~ nig /nip in each pair-producing colUsion we can 
treat the losses as a continuous process and write 



Lp"" = -g:^[ypnp{yp,myp,t)\ 



(22) 
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where 



(23) 



is the convolution of the normahsed colUsion rate with the in- 
elasticity ^ given by eqn. (19). 



4. Proton injection and blacl(-body photon field 

As a first example, we examine the case of relativistic proton 
injection and subsequent cooling on a black-body photon field. 
We assume, as usual in astrophysical cases, a power-law proton 
injection of the form 



2p(r) = Gp.or ^^(r - rp,min)^/(rp,max - r) 



(24) 



where yp^min and yp,max are respectively the lower and upper cut- 
off of the proton distribution. We also assume that the protons 
can escape at a rate fp '^.^ from the source region, taken to be a 
sphere of radius R. Then the stationary solution, in the absence 
of proton losses is simply 



«p(r) = Gp.O^p.escT ^Hiy - rp,min)^(ri 



p.max 



r). 



(25) 



We assume next that a significant number of external pho- 
tons is present and that these have a black-body distribution 
of temperature Tbb- The energy density Ubb of these photons 
in the source can be parameterised in terms of the black-body 
compactness /'bb, defined by the relation 



^BB = 



UbbO'tR 



This emission can arise, for example, from the surface of an ac- 
cretion disk located close to the source. In general, if the source 
is irradiated by a black body at temperature Tbb whose surface 
occupies a soUd angle of A£l as seen from the source, then 



Ubb = l^lflradrBB 



(27) 



where Onid is the radiation density constant. In terms of com- 
pactness; 



^BB = 615 X 



AQ 
4n 



jXlO^Kl VlOi^cm^ 



(28) 



We further assume that initially there are no other pho- 
tons or electrons present (except maybe in very small num- 
bers) and that there is no external electron injection. Since 
the Bethe-Heitler pair production has a threshold condition 
jx > I + nif^/nip (corresponding to head-on collisions) it fol- 
lows that, in the case where y = 7p,max(^rBB/weC^) «: 1, there 
will be negUgible pair production. Therefore, for all practical 
purposes, the solution of the proton equation will still be given 
by Eq. (25). However, for J ~ 1, the proton-photon reaction rate 
is substantial and cannot be neglected as a proton loss mecha- 
nism. Moreover, the produced electron/positron pairs provide 
an injection term for the electron equation and lose energy 
mainly by inverse Compton scattering and/or synchrotron ra- 
diation. As a result a photon spectrum is formed. To investigate 
the properties of this emission, such as its luminosity and spec- 
tral shape etc., we must distinguish between two cases: 




logx, log7^ 

Fig. 3. Steady-state photon spectrum resulting from a power-law pro- 
ton injection and subsequent proton-photon pair production on a 
black-body photon field in the case where photon-photon pair pro- 
duction (i) has been ignored (dashed line) and (ii) has been taken into 
account (full line). The dotted line curve depicts the electron/positron 
distribution function at production. For this particular run the protons 
were assumed to be injected with a power-law of slope j8 = 2, between 
the Umits yp.min = 10° ' and yp.max = 10*. Also the values {p = 2.5 10'^ 
(26) and fp.csc = have been assumed. The black-body photon field pa- 



rameters were Tbb = 10^ K (0 : 
losses were neglected. 



1.7 10"^) and 4b = 1- Synchrotron 



4.1. Case 1:{p ^ ^bb 3nd negligible synchrotron 
losses 

In this case, the photons produced by electrons created in the 
BH process are not important as targets and the resulting pho- 
ton spectrum is quite simple. It is a single power-law of (num- 
ber) spectral index -3/2, extending up to an energy Xmax, that is 
approximately the inverse of the temperature of the black-body 
field, i.e. Xn,ax - © where = (kTsB/meC^)- The explana- 
tion of this spectrum is straight-forward: protons pair-produce 
on the black-body field and, since they are produced with 
high energies, the pairs cool on the black-body photons ini- 
tially by Compton scattering in the Klein-Nishina regime. This 
produces y-rays that are above the threshold for pair produc- 
tion on the black-body photons and, therefore, are re-absorbed. 
Lower energy pairs, which cool by Compton scattering in the 
Thomson regime, produce photons which are below the thresh- 
old for photon-photon pair production. This naturally produces 
an electron distribution function oc y"^ and thus a photon 
spectrum Uy oc x~^^^. The value of Xmax. therefore, is set by 
the condition Tyy(xn,ax) - 1^ where Tyy is the optical depth for 
photon-photon pair production on the background black-body 
field. 

Since for the present case there is neither escape of elec- 
trons from the system nor any sink of energy, (e.g. synchrotron 
self-absorption) other than photon escape, the radiated lumi- 
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nosity (or, equivalently, equals the luminosity injected in 
pairs once a state steady is achieved, i.e. ty ^ where 



(29) 



is the compactness of the created pairs from the Bethe-Heitler 
pair production. Therefore, the overall photon luminosity de- 
pends on the parameter y defined above and, as we show below, 
also on ^BB and t^. 

As an illustrative case we take ^bb = 1 » = 2.5 10"^. 
Here also Tbb = 10^ K, while rp,„u„ = 100 \ yp.max = 10^ 
The dotted line in Fig. 3 depicts the electron injection func- 
tion which shows a broad maximum. The same figure shows 
the photon spectrum which is obtained from the cooling of 
these electrons and the ambient black-body field. The dashed 
line shows the spectrum in the case where the yy pair produc- 
tion has been artificially switched off. The spectrum is flat and 
peaks at high energies. We note that both the electron produc- 
tion function and the unabsorbed photon spectrum extend to 
about two orders of magnitude above jp^mx-, an effect that is due 
to the kinematics of the Bethe-Heitler pair production - see Fig. 
3. At lower energies it produces the characteristic -3/2 power- 
law. The full line shows the final photon spectrum which in- 
cludes yy absorption. It has still the same power-law, however 
all the details of the electron injection spectrum which were 
evident in the unabsorbed spectrum have disappeared due to 
the intense attenuation. Nevertheless the overall luminosity is 
conserved (as it should be) and as the electromagnetic cascade 
redistributes the power to lower energies, the flux is increased 
there. Figure 4 shows the effects of losses on the steady-state 
proton distribution for these parameters (/"bb - 1) and two 
higher values of ^bb, compared to the loss-free case. 

The efficiency of the BH process, i.e., the ratio of the proton 
power turned into pairs to the total power injected as protons, 
is shown in Fig. 5 as a function of the maximum Lorentz factor 
of injection 7p,max- The inelasticity of this process is small (f ~ 
10"^), so that high efficiency can only be achieved if a proton 
interacts many times before escape. This is indeed the case for 
high black-body compactnesses — the efficiency exceeds 40% 
for /"bb > 10"* and yp.max > 10^. 

4.2. Case 2; <^ ^bb ^nd non-negligible synchrotron 
losses 

In the more realistic case, where synchrotron losses cannot 
be neglected, the spectrum of the electrons becomes more 
comphcated, as now the electrons cool by a combination of 
synchrotron radiation and inverse Compton scattering. Fig. 6 
shows the obtained spectra in the cases where the magnetic 
compactness ^b> defined according to 



{b = cnR 



(30) 



is comparable to ^bb- The redistribution of the radiated lu- 
minosity to increasingly lower energies is evident, the rea- 
son being that, for this particular choice of T^b and B, in- 
verse Compton scattering produces much harder photons than 




Fig. 4. Proton .steady-state spectra for various external black-body 
photon compactnesses. All cases are taken for the same injection pro- 
ton parameters (^p = 2.5 10"^, yp,max = 10*, fp,esc = 1) and for the 
same temperature of the external black-body field (rBB = 10' K) The 
plotted proton spectra are shown when = (no-loss case, full line), 
4b = 1 (dotted line), 4b = 10^ (short-dashed Une) and 4b = 10^ 
(long-dashed line). 



the synchrotron mechanism. Therefore, as €b progressively in- 
creases, the peak is shifted towards the lower energies, showing 
the increasing importance of synchrotron radiation as an elec- 
tron energy loss/photon emission mechanism. As in Fig. 3, the 
spectrum continues to be strongly absorbed for photon energies 

X > Xmax — ® 



5. The Pair-Production/Synchrotron Instability 

5.1. Marginal stability 

Kirk & Mastichiadis (1992, henceforth KM92) have shown 
that ultra-relativistic protons can, under certain conditions, be- 
come unstable to various types of radiative instability. They 
showed expUcitly the necessary conditions for one of them 
to occur, namely the Pair-Production/Synchrotron instability 
(henceforth PPS). To understand the basic idea, assume that 
protons are confined in a region of characteristic radius R 
where a magnetic field of strength B is also present. Assume, 
moreover, that the protons are relativistic and have Lorentz fac- 
tors such that if they photo-pair produce, the synchrotron pho- 
tons radiated from the created pairs are sufficiently energetic 
for the protons to produce more pairs on them. Making the sim- 
pUfying assumptions that (i) the created pairs have the same 
Lorentz factors as the protons and (ii) the synchrotron photons 
are all emitted at the critical frequency, KM92 showed that in 
order for protons to be able to initiate this loop they should have 
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Fig. 5. Proton efficiency, i.e. fraction of power lost to pairs to total 
power injected in protons as a function of the upper cutoff of the proton 
distribution -yp max- The rest of the proton injection parameters are kept 
constant (f,, = 2.5 10"', 7p.n,i„ = 10"', fp^sc = 1), while the external 
black-body field has temperature Tbb = 10^ K. and compactnesses 
(bottom to top) 4b = -01, 1, 100, 10". 

Lorentz factors above a critical value given by 

rcri. = (-] (31) 

where b - BIB„ with the critical value of the magnetic 
field (Bcr = mlc^/eh = 4.413 lO^^ Gauss). For this loop to be 
self-sustained it is necessary that at least one of the synchrotron 
photons should produce a pair before escaping the source and 
this condition naturally leads to a critical proton number den- 
sity (see KM92, eq. 6). 

MK95 presented a numerical simulation of the PPS insta- 
biUty in the case where protons are accelerated from low mo- 
menta by a Fermi-type acceleration scheme. Once the protons 
(assumed to have a density exceeding the critical number den- 
sity) reached energies above ycrit> the conditions for the insta- 
bility loop were complete and the internally produced photons 
increased, saturating the acceleration, and driving the system to 
equilibrium. However, the code used by Mastichiadis & Kirk 
(1995) is limited by the simplifying assumptions mentioned 
in the previous paragraph. Here we re-examine this problem 
with the improved version of the code that uses, as described 
in section 3, the Bethe-Heitler pair-production spectra as given 
by Monte-Carlo code and the full synchrotron emissivity. The 
objective is to find an accurate estimate of the critical proton 
number density above which the PPS instability occurs, i.e., a 
numerical version of the analytical (but approximate) Fig. 1 of 
KM92. 

To make our results directly comparable with those of the 
aforementioned figure, we have set the parameters to the val- 
ues prescribed there. For this we took a source size R = 10^^ 




logx 

Fig. 6. Steady state photon spectra for the same parameters used in 
Fig. 3 but including synchrotron radiation. The magnetic field used 
was B = 100 G and the magnetic compactnesses were £b = 0.01 (full 
line), /"b = 0.1 (dashed Une) and ^b = 1- The black-body compactness 

was £bb = 1- 

cm, a magnetic field B = 10^ Gauss and a proton distribution 
function of the form ripijp) = n^fij^^ with (3-2 from the low- 
est allowed proton energy 7p,min = 10*^ ' to a maximum energy 
Tp.max- We note that in this case the proton distribution is held 
constant throughout each run, i.e. protons do not evolve. For 
various values of yp.max, we run the code for different values of 
the only remaining free parameter (np,o). According to KM92, 
the time evolution of the photon and electron distribution func- 
tions is of the form ny{t) oc ri^ e^' with s > when the protons 
are in the unstable regime. Thus in order to verify numerically 
the existence of the PPS instability we seek a value n™,' above 
which the internally produced electron/positron pairs and pho- 
tons start increasing with time. 

The onset of instability for jp^max = 10'' can be seen in 
Fig. 7 which depicts the photon compactness ly as a function of 
time t (expressed, as always, in units of ta) for various values 
of Hpfl around «p"'. When the protons are in the stable regime 
there is some pair production between the protons and the syn- 
chrotron photons produced from the initial electron distribution 
but the system, for times larger than the synchrotron cooling 
time, settles to a steady-state. Thus for t > /cooi we get s = 0. 
However, as can be seen from the figure, as Wp o increases (from 
bottom to top) the photons start to grow exponentially, with s 
increasing with increasing ripfi. It is worth mentioning that each 
curve corresponds to a value of np,o that is larger by only 2% 
than its previous value; therefore, this figure depicts the rapid 
onset of the instability. 

Fig. 8 shows the behaviour of * as a function of the proton 
normalisation Wp o for three values of 7p,max- It is evident that 
once the instability sets in, s is a very sharp function of «p_o. 
Therefore, we find that increasing approximately by a fac- 
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Fig. 7. Plot of the internally produced photon compactness £y as a 
function of time for values of Mp,o around the critical value n^Q. In 
each curve the value of ti^o is increased 2% from its previous value. 
For these runs yp,max = 10*. 
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Fig. 8. Behaviour of the instability growth index s as a function of n^p 
for values of 7p,max = 10^ 10^ 10' (right to left). 

tor of 2 above its critical value, s becomes greater than one, 
i.e. the density of photons starts growing on a timescale shorter 
than the crossing time of the source. This, as we shall see in the 
next section, has catastrophic consequences for the high energy 
protons as the spontaneously growing photons make them lose 
their energies. 

Fig. 9 shows the marginal stability curve as is obtained from 
the present code (in practice we have calculated the values of 




4 5 6 



log(7p.„,ax) 

Fig. 9. Plot of the numerically obtained marginal stability as a function 
of 7p,max (full line). The dashed line shows the i = 1 locus, while the 
dotted line depicts the marginal stability results of KM92. Note that 
the KM92 approximation is reasonably good at low yp.max- 



np,o which correspond to s = .05. Due to the steep dependence 
of s on np_o this value can be considered as very close to the 
marginal stability one). Note that this curve is in very good 
agreement with the curve estimated by KM92 (plotted here as 
a dotted Une) for values of y close to the threshold, but exceeds 
it by a factor of about 2 for at high energy. This difference can 
be attributed to the overestimation of the electron production 
rate in the BH process as discussed in Section 2. For maxi- 
mum values of the proton Lorentz factor close to ycrit, all the 
BH interactions occur close to threshold, so that the assump- 
tion used by MK is accurate. (The small difference in the shape 
of the KM curve and our present result in this energy range 
can be attributed partly to the kinematics of Bethe-Heitler pair 
production and partly to the fact that the full expression for 
the synchrotron emissivity was used). However, once the upper 
cut-off of the proton distribution substantially exceeds y„ii, the 
effect of pairs injected with y > becomes important. These 
are not taken into account in the approximation used by MK92, 
but are treated accurately in the simulation-based method used 
here. It is evident that the basic concept of the feedback loop 
remains unaltered by our more accurate treatment. The quanti- 
tative impUcations, however, are discussed below. 

5.2. Proton injection 

In order to see the effects of the PPS instability when proton 
losses are taken into account we assume once again (see section 
4) that protons are injected in a region of radius R with a power- 
law (eq. 24). However, to make the picture less complicated, 
we assume that there is no external photon field present. Thus 
one expects that the proton spectrum will reach an equilibrium 
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state given by the no-loss solution (eq. 25). This is indeed trae 
as long as this steady-state solution is below the critical density 
for the PPS instability, i.e. when Wp^o = gp.ofp.esc < np"J(rp,max), 
with ripfi the normalisation of the protons and fip'QCyp.max) given 
from Fig. 9. In this case also there is a very small number of 
internally produced photons. However, as soon as the condi- 
tion tipfi = 2p,ofp,esc > «p"o(rp,max) holds, the Criteria for the 
PPS instabiUty are satisfied. As a result, the photons increase 
exponentially, protons lose energy due to pair production and 
the system behaviour depends on the choice of the parameters 
Qpfi and tp^esc- 

As a first example, we show in Fig. 10 the case where 
fi - 2, yp,max = 10*, ?p,esc = 1 with a proton injection rate 
of 2p,o = 0.22. The proton compactness parameter, according 
to eq. 9, is ^p - 1710. The above combination of 2p_o and fp,esc, 
according to Fig. 8, corresponds to a feedback loop that causes 
the photon density to increase with s = 0.7. The dotted and 
short-dashed lines show, respectively, the evolution of the low- 
est and highest differential density bins of the proton energy. 
In agreement with the analytical solution of Eq. (5) in the no- 
loss case, the number of particles in these bins increases very 
quickly (in about one 1^) to a steady state, which, however, cor- 
responds to an unstable proton configuration. The long-dashed 
line shows the photon compactness when only synchrotron ra- 
diation is taken into account. This increases as {y oc e^' with 
s = 0.7 until it reaches a steady state. The dot-dashed line 
shows the effects on the photon compactness caused from the 
addition of inverse Compton scattering and photon-photon pair 
production. As these processes are quadratic with respect to the 
photon and electron/positron number densities they do not af- 
fect the slope of the compactness at early stages, i.e. as long 
as <sc 1. However, at the later stages of evolution these 
processes become important and, because they are quadratic, 
they cause the photon compactness to increase even faster and 
reach saturation sooner. Finally the solid line shows the evolu- 
tion of the photon compactness when, in addition to the above 
processes, photon trapping due to the high density of created 
pairs is taken into account. This leads to higher photon com- 
pactnesses as relatively more energy is extracted from the high 
energy protons. This curve must be considered the 'correct' one 
as it contains all the relevant processes. On the other hand the 
photon increase causes the high energy protons (short-dashed 
line) to lose energy and settle in a new steady state. These 
losses do not affect naturally the low energy protons (dotted 
line) which maintain their original steady state. 

At the other extreme, one can envisage a case where pro- 
tons are injected slowly, but have a very long escape time. In 
this case we take 2p,o = .22 10"^ and ?p,esc = 10^ which corre- 
sponds to the same no-loss steady state as before. The result is 
shown in Fig. 11, where quasi-periodic behaviour typical of a 
relaxation oscillator is apparent. The protons accumulate in the 
source and, once their density rises above the critical density, 
the photon density grows rapidly (on the timescale of a few 
times tcr) and deplete the high energy part of the proton spec- 
trum. Once the protons have lost their energy, there is nothing 
left to sustain the loop and the photons escape from the system. 
This cycle is repeated as the protons accumulate again in the 



4 




' 1 


1 1 1 1 1 1 1 


2 




i // 

if 1 


- 






1 1 








1 








1 









1 


- 




'' 


J \ 
f ' 




2 




/ 1 

/ \ 

1 \ 
/ \ 




4 








6 




1 





50 100 
t 

Fig. 10. Plot of the evolution of the system in the case where protons 
are injected with p = 2, yp,max = 10*, /p.esc = 1> while the combina- 
tion of the proton injection rate (gp,o = 0.22) and fp esc is such that the 
steady-state proton distribution in the no-loss case corresponds to an 
unstable proton distribution. The solid line shows the photon compact- 
ness when all processes are taken into account, while the long dashed 
and dot-dashed lines show the photon compactness when certain pro- 
cesses are omitted (for details see text). The dotted and short dashed 
lines show the evolution of the first and last proton occupation number 
bin in the case where all the relevant processes are taken into account. 



source. The behaviour of the lowest and highest proton energy 
bins is also shown (dotted and dashed lines respectively). It is 
clear that photons and high energy protons are anticorrelated, 
in the sense that when one population is high, the other is low. 
These cycles are similar to those found by Stem & Svensson 
(1991) using Monte-Carlo techniques. 

6. Simultaneous electron and proton injection 

As a last case we examine the situation where electrons and 
protons are injected simultaneously, i.e., both proton and elec- 
tron kinetic equations (5) and (6) have an external injection 
term. We assume that both of these terms are in power-law 
form, so that, in addition to Eqn. (24) that describes proton 
injection, we prescribe electron injection using a similar ex- 
pression: 

(re) = Q.,07'^H(y, - re,mm)//(re,max " Te)- • (32) 

To lower the number of free parameters we assume that 
both electrons and protons are injected with the same spectral 
index p and with the same maximum Lorentz factor, i.e. we set 
arbitrarily 7e,max = 7p,max- Moreover, when solving the relevant 
kinetic equations, we assume that the two species have equal 
escape times, i.e., fp esc = ^e.esc- 

Depending on the particular choice of parameters, high en- 
ergy electron injection can result in a synchrotron and/or an 
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Fig. 11. Plot of the evolution of the system in the case where a low 
proton injection rate (2p,o = 2.2 10 combines with a slow proton 
escape time (?p_esc = 10^) in such a way that the steady-state proton 
distribution in the no-loss case to correspond to an unstable proton 
distribution. The rest of the parameters are as in the previous figure. 
Solid line shows {y,, while dotted and short dashed lines show the evo- 
lution of the proton spectrum at 7p = 10." ' and 7p = 10'' respectively. 



inverse Compton component (see, for example, Mastichiadis 
& Kirk 1997). When protons are injected as well, these will 
interact with the aforementioned photons causing a secondary 
injection of Bethe-Heitler pairs. This leads to a non-linear sit- 
uation. To see this, one should compare the compactness of 
the externally injected electrons £f'^ (eq. 10) with the corre- 
sponding compactness of the internally produced pairs via the 
Bethe-Heitler process (Eq. 29). As the latter is, in general, 
a function of both {p and ^g"', when < ^™ the system op- 
erates in the non-linear regime, in the sense that the cooling of 
the protons occurs mainly on the internally produced photons. 

To investigate the effect described above we proceed as fol- 
lows: we keep the electron injection parameters constant and 
change only the normalisation of the proton injection rate Qp o 
in eq (24) — or, equivalently, we change the proton compact- 
ness parameter £p (Eq. 9). When (p <k £f- the resulting pho- 
ton spectrum is simply that produced from the cooling of the 
externally injected electrons. However as we increase ip, the 
quantity €™ increases as well and, at some stage, it becomes 
comparable to {f^K Above this point the system enters the non- 
hnear regime, as the cooling of the protons occurs primarily 
on its own radiation. Finally, above some critical proton com- 
pactness, a loop analogous to the one described in Section 5 
operates, and the protons convert a substantial fraction of their 
energy content to electron/positron pairs and radiation. This 
is depicted in Fig. 12 which shows the evolution of the pho- 
tons when ^^"^ is kept constant, while {p varies. The first curve 
(dotted Une) assumes that no protons are injected. Due to the 



particular choice of the parameters the primary electrons cool 
fast and the system quickly reaches a steady state. In the next 
two cases, the injected proton compactness is ip = 200 (long- 
dashed line) and (p = 400 (short-dashed line) respectively. 
Again a steady state is quickly reached. However the photon 
compactness {y increases as {p increases because of the BH 
pairs created and subsequently cooled. In both of these cases 
the feedback loop does not operate in the sense that the proton 
losses remain low — or, equivalently, iy ^ ^g"' + i^^ <s; {p. 
However, for the last two cases which are for (p = 800 (full 
line) and dp - 1600 (dot-dash line) the effects of the feedback 
are evident. In the case of the solid curve, and for the first 300 
or so Ught crossing times, there is a gradual increase of pairs 
and photons in the system until their numbers are built to a 
level that allows a catastrophic release of the energy stored in 
protons. After that, a steady state is reached, but at a level much 
higher than that of the previous cases: fy ^ /'g"' + ^™ ^ {p. It is 
worth mentioning that this loop operates at a proton compact- 
ness which is below the critical threshold obtained in Section 5. 
Therefore, the presence of external electrons helps to initiate 
the catastrophic proton energy losses at lower proton densities. 
Finally, the uppermost curve corresponds to a proton injection 
which is above the critical threshold for the PPS instability and 
the photons grow very quickly as discussed in Section 5 — see 
also Fig. 10. 



The photon spectra corresponding to the steady states ob- 
tained in each of these runs are shown in Fig. 13. The bot- 
tom curve (dotted hne) corresponds to the injection of elec- 
trons only whereas the four others are for electron and pro- 
ton injection, corresponding to the curves of Fig. 12. The ex- 
tra component due to BH pair production is apparent on the 
two lower spectra that include protons (long and short-dashed 
lines). The two uppermost curves correspond to the steady-state 
spectra when the loop was able to extract most of the proton en- 
ergy. Due to the resulting high photon compactness the spectra 
are strongly absorbed at energies above 1 MeV due to photon- 
photon pair production. 



The above findings are summarised in Fig 14 which shows 
the photon compactness of the system versus the proton com- 
pactness for various injected electron compactnesses. For low 
values of £p (i.e. less than 30), the photons of the low-frequency 
part of the SED produced come almost exclusively from the 
presence of the electrons, i.e. the protons cannot significantly 
influence the low-frequency behaviour of the system. We note 
here that for synchrotron proton blazar model fits to the SED of 
BL Lac Objects observed at gamma-ray energies (e.g. (Miicke 
et al. 2003)) the proton compactness is {p ~ 10^^-10"^ im- 
plying that Bethe-Heitler pair production, and any associated 
instability, is unimportant in these models. For intermediate 
values of {p (i.e. between 30 and 300) the internally produced 
Bethe-Heitler pairs make their presence visible in that their 
cooling increasingly dominates the photon spectrum. Finally 
at even higher compactnesses the protons become unstable and 
cool efficiently on their "own" radiation. 
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Fig. 12. Plot of the evolution of the photons in the case where protons 
and electrons are injected simultaneously with the same slope p = 2 
and the same high Lorentz factor cutoff, i.e. yp.max = Tcmax = 10* 
in a region of radius R = 10'^ cm immersed in a magnetic field of 
strength B = 10^ G. The injected electron compactness is = 0.04 
and the evolution of the photons is shown for proton compactnesses 
(bottom to top) = 0, 200, 400, 800 and 1600. Both species have 
an escape time equal to AH except the highest curve correspond to 
stable proton distributions in the no-loss case. Nevertheless the solid 
curve shows a loop that eventually leads to catastrophic proton losses. 
The uppermost curve (dot-dashed line) corresponds to a proton distri- 
bution that is unstable even in the no-loss case. Therefore the photons 
increase quicldy and drive the system to an equilibriimi. 



7. Summary/Discussion 

In the present paper we have examined some consequences 
arising from the presence of ultrarelativistic hadrons in com- 
pact sources. This was done with the help of a numerical code 
that was constructed to follow the evolution of the system 
through the solution of three coupled, time-dependent kinetic 
equations for protons, electrons and photons respectively. All 
the relevant basic processes involving electrons and photons 
in astrophysical pair-plasmas were included. The coupling be- 
tween the hadronic and the leptonic component was assumed to 
occur via Bethe-Heitler pair-production. For this process, de- 
tailed electron/positron pair-production spectra were obtained 
with the help of a Monte-Carlo code. These were then incor- 
porated into the kinetic equations which were subsequently 
solved numerically, revealing effects mainly due to synchrotron 
and inverse Compton losses. 

The choice of the kinetic equation approach allowed us to 
study various aspects of the behaviour of such a system. Thus 
we showed that the presence of an external black-body radi- 
ation field can extract energy efficiently from the relativistic 
protons only when the photon compactness is high. 

Of special interest are the non-linear cases, i.e., cases in 
which the protons cool not on a prescribed external photon 



Fig. 13. Plot of the photon spectra corresponding to the steady states 
obtained for the runs shown in Fig. 12. The bottom curve (dotted) 
shows the spectrum obtained from the cooling of the injected elec- 
trons when no protons are injected. As the injected proton compact- 
ness increases the photon spectrum is modified due to the presence of 
the radiation from the created pairs (long and short-dashed curves). 
Finally when the loop operates the photon spectrum is strongly ab- 
sorbed above 1 MeV by photon-photon pair creation (solid and dot- 
dashed curves). 




Fig. 14. Plot of the photon compactness versus the proton compactness 
for various injected electron compactnesses. The curves correspond to 
= 1, 10 and 100 (from bottom to top). The other parameters are the 
same as those used in Figs. 12 and 13. The dotted line represents the 
critical proton compactness in the no-loss case as this was described 
in Section 5. 
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field, but on the radiation of the internally produced Bethe- 
Heitler pairs. In the present paper we have verified the ex- 
istence of the 'Pair-Production/Synchrotron' loop previously 
studied analytically. We showed that this is a very eflftcient way 
of channelling proton energy into electron/positron pairs and 
radiation, and that the coexistence of relativistic electrons in 
the system does not stabilise the system but, on the contrary, 
lowers the critical density threshold, i.e. it facilitates the effi- 
cient transfer of energy from the hadronic component to the 
leptonic/photonic one. This, and the fact that the threshold and 
critical density conditions can be greatly relaxed if the protons 
are in relativistic bulk motion, makes this loop a promising 
candidate for some AGNs (Kazanas & Mastichiadis 1999) and 
GRBs (Kazanas et al. 2002). 

Of the various processes omitted in the present treatment, 
the most important is that of photo pion-production. This, how- 
ever, does not aff'ect the results of the present paper because, in 
the examples shown, the initial conditions were chosen so as to 
avoid the onset of this process. Other processes involving pro- 
tons, such as proton synchrotron radiation and proton-proton 
interactions, etc., are negligible for the parameters of the par- 
ticular examples given in this paper. 
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